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THE DYNAMICS OF RADIATIVE SHOCK WAVES: LINEAR AND 

NONLINEAR EVOLUTION 

A. Mignone 1 ' 2,3 
ABSTRACT 

The stability properties of one-dimensional radiative shocks with a power-law cooling function 
of the form A oc p 2 T a are the main subject of this work. The linear analysis originally presented 
by Chevalier & Imamura, is thoroughfully reviewed for several values of the cooling index a and 
higher overtone modes. Consistently with previous results, it is shown that the spectrum of the 
linear operator consists in a series of modes with increasing oscillation frequency. For each mode 
a critical value of the cooling index, a c , can be defined so that modes with a < a c are unstable, 
while modes with a > a c are stable. The perturbative analysis is complemented by several 
numerical simulations to follow the time-dependent evolution of the system for different values 
of a. Particular attention is given to the comparison between numerical and analytical results 
(during the early phases of the evolution) and to the role played by different boundary conditions. 
It is shown that an appropriate treatment of the lower boundary yields results that closely follow 
the predicted linear behavior. During the nonlinear regime, the shock oscillations saturate at a 
finite amplitude and tend to a quasi-periodic cycle. The modes of oscillations during this phase 
do not necessarily coincide with those predicted by linear theory, but may be accounted for by 
mode-mode coupling. 

Subject headings: hydrodynamics - instabilities - methods: numerical - shock waves - stars: 
binaries: close 



1. INTRODUCTION 

Radiative shock waves are believed to play a key role in a variety of different astrophysical environments, 
including magnetic cataclysmic variables (Wu 2000; Cropper 1990), jets from young stellar objects (Hartigan 
et al. 1994), magnetospheric accretion in T-Tauri stars (Calvet & Gullbring 1998), colliding stellar winds 
(Stevens et al. 1992; Antokhin et al. 2004) and supernova remnants (Kimoto & Chernoff 1997; Blondin et 
al. 1998; Walder & Folini 1998). 

Most of the earlier theoretical investigation has been motivated by the dynamics of accreting shocks in 
magnetic cataclysmic variables. In these systems, a strongly magnetized white dwarf (10 6 < B < 10 8 Gauss) 
accretes matter directly from a late-type star without the formation of a disk. Instead, the mass transfer 
process is magnetically channeled and matter accumulates through a stand-off shock on a small fraction of 
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the stellar surface. At temperatures of 10 8 -j- 10 9 K, the shock-heated plasma radiates its energy via optically 
thin bremsstrahlung becoming a powerful source of X-rays. 

It has been shown that, under certain circumstances, the postshock flow is subject to a global thermal 
instability (more precisely, an overstability) caused by rapid variations of the cooling time scale with the 
shock speed. The instability drives the shock front to oscillate with respect to its stationary position, causing 
variations in the amount of radiation emitted from the postshock region. The instability mechanism has 
been invoked in the past to explain the optical quasi-periodic-oscillations (QPO) observed in AM Her-type 
systems (Larsson 1992; Middleditch et al. 1997; Wu 2000). Similarly, a relevant issue arises in questioning 
the validity of steady shock models with shock velocities v s > 130 km s _1 , routinely used in interpreting 
emission line observations from interstellar shocks (Innes et al. 1987a, b; Gaetz et al. 1988; Sutherland ct al. 
2003). 

The nature of the instability was first studied analytically by Chevalier & Imamura (1982, CI hereafter), 
who presented a linear stability analysis of planar radiative shocks with volumetric cooling rate A oc p 2 T a . 
CI showed that the shock has multiple modes of oscillation, and the stability of a particular mode depends 
on the value of the cooling index a. In general, higher power dependences on the temperature were shown 
to inhibit the growth of instability. Thus the fundamental mode was shown to become unstable for a < 0.4, 
while the first and second harmonics are de-stabilized when a < 0.8. Higher order harmonics were not 
considered by CI. 

Perturbative studies of one-temperature flows with a power-law cooling function were afterwards con- 
sidered by several authors. Imamura (1985) presented a linear stability study of radiative shocks where the 
cooling function had a weaker dependence on density, i.e. A cx pT a . Bcrtschinger (1986) examined the 
structure of spherical radiative shocks and also considered the effects of nonradial perturbations. He showed 
that modes that are stable to radial perturbations may become unstable for small transverse wavenumbers. 
For sufficiently large wavenumbers, however, all modes are eventually stabilized. 

Effects due to cyclotron emission were considered in Imamura et al. (1991), Wu et al. (1992), Wu et al. 
(1996). Noise driven models were proposed by Wood et al. (1992). Effects of gravity were studied by Houck 
& Chevalier (1992), while magnetic field effects were considered by Edelman (1989a), Edelman (1989b), Toth 
& Draine (1993), and by Hujeirat & Papaloizou (1998) in the two-dimensional case. Dgani & Soker (1994) 
considered radiative shocks with a mass loss term. 

Recently, Yamada & Nishi (2001) investigated the stability properties of shock-compressed gas slabs by 
introducing a cold layer of finite thickness. They considered both symmetric and antisymmetric modes and 
reported the existence of quasi-oscillatory modes, in addition to the overstable ones. 

Stability properties of radiative shocks with unequal ion and electron temperatures were investigated 
by Imamura et al. (1996), while flows with multiple cooling functions (i.e. bremsstrahlung and cyclotron) 
were examined both in the single and two-temperature regimes by Saxton et al. (1997), Saxton et al. (1998), 
Saxton (1999), Saxton & Wu (1999), Saxton & Wu (2001) and Saxton (2002). Additional references can be 
found in the review by Wu (2000). 

Investigation of the full time-dependent problem has also received extensive attention over the past two 
decades and several numerical simulations have been carried along. The oscillatory instability was discovered 
in the first place numerically by Langcr ct al. (1981, LCSa hereafter) and Langer et al. (1982, LCSb hereafter), 
who investigated spherically symmetric accretion onto magnetized white dwarfs. They showed that flows 
with volumetric cooling rates A ~ p 2 T a become unstable when a < 1.6, a limit subsequently revised to 
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a < 0.6 in Langer et al. (1983). In the case of bremsstrahlung cooling (a = 1/2), the shock position was 
shown to undergo periodic oscillations over the surface of the white dwarf. Imamura et al. (1984, IWD 
hereafter) and Imamura (1985) used a Lagrangian code to investigate the stability of radiative shocks with 
a power-law cooling function A ~ p 2 T a . Effects of gravity, thermal conduction and unequal electron and 
ion temperatures were also included. In IWD, the critical values of a (above which oscillations are damped) 
were shown to lie somewhere in the range 1/3 < a < 1/2 and 1/2 < a < 0.6 for the fundamental and 
first-overtone, respectively. Wolff et al. (1989) considered both one and two-temperature calculations; the 
one-temperature calculations showed that the fundamental mode is unstable for a = 0.33, first and second 
harmonic are unstable when a = 0.65, while the system is stable at a = 1. The one-dimensional calculations 
of Strickland & Blondin (1995, SB hereafter) showed that for flows incident into a wall, large amplitude 
oscillations are damped when a > 0.5. SB also considered perturbed steady-state models, showing that 
power-law cooling functions with a < 0.75 produced shock oscillations which saturated at a finite amplitude. 
The results of SB were recently extended by Sutherland et al. (2003, SBD hereafter) to a more realistic 
cooling function. 

Most results of the previous numerical investigations, however, bear no clear relation with the predicted 
linear behavior, and a direct comparison with perturbative studies has proven to be only partially successful. 
Although the salient features of thermally unstable shocks have been commonly reproduced, the modes 
of instabilities can not always be identified with the linear ones, with the exception of one or two modes. 
Besides, controversies exist on the value of the critical a above which the system should become stable. 
These apparent inconsistencies may be partially due to the fact that most calculations do not include a 
stationary solution in their initial condition, which makes the problem inherently nonlinear since the very 
beginning. Moreover, a substantial disagreement exists between Eulerian and Lagrangian calculations and 
on the numerical treatment of the lower boundary condition which plays a crucial role in the dynamics of 
the post shock flow. 

Here I wish to present some new and detailed calculations in the attempt to settle part of these con- 
troversies. The results of this work will also serve as a basis for future extensions, where effects of magnetic 
fields and more realistic cooling functions will be considered. 

The paper is organized as follows. In §2 the problem is defined and the relevant equations are introduced. 
In §3 the stability properties of 1-D planar radiative shocks are reviewed in more detail, while in §4 numerical 
simulations are presented with particular emphasis to the comparison with linear theory and to the choice of 
boundary conditions. A new time-dependent treatment of the lower boundary is introduced and the details 
of implementation are given in Appendix A. 

2. STATEMENT OF THE PROBLEM AND EQUATIONS 

Consider a one-dimensional supersonic flow with uniform density p m and velocity V[ n , initially propa- 
gating in the negative x-direction, i.e., V{ n = — \vi n \. The flow is brought to rest by the presence of a rigid 
wall located at x — 0, and a shock wave forms at some finite distance x s from the wall. Through the shock, 
the bulk kinetic energy of the incoming gas is converted into thermal motion and the flow decelerates to 
subsonic velocities. In the postshock region, the thermal energy of the accreting gas is radiated away by 
cooling processes. 

In steady-state, the dynamical time scale is equal to the cooling time scale, so a fluid element travels 
through the postshock region and cools exactly to zero temperature by the time it reaches the wall. 
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In several areas of interest and to make the problem more tractable, radiative losses are treated in the 
optically thin regime and the cooling rates are normally specified as functions of the temperature, density 
and relative abundances. In this work it is assumed that the volumetric cooling rates can be described by a 
single power-law (in temperature) function 

A(p,p) = ap 2 (^y , (1) 

where a is a physical constant depending on the particular cooling process, a is the cooling index and p and 
p are, respectively, the density and pressure of the gas. 

The problem of a supersonic flow into a wall is, of course, a simplified abstraction of a more complex and 
specific physical setting. Effects due to thermal conduction, magnetic fields and multi-dimensional effects 
are neglected in this paper and will be considered in future works. With these assumptions, the problem 
can be described by the Euler equations for a one-temperature flow in planar geometry: 

dp dv dp 

i +p d- x +v £ = ^ (2) 

dv ^ dv ^ 1 dp q 
dt dx pdx 

where v is the fluid velocity, C is a constant and an ideal equation of state with constant specific heat ratio 
7 has been assumed. Equations (2) through (4) are put in a dimcnsionless form by expressing density and 
velocity in units of their inflow values, i.e., p m and \v ia \. With this choice, the flow variables immediately 
ahead of the shock become p = 1, v = — 1 and p = l/{^fM 2 ) 1 with M being the upstream Mach number. 

The length scale of the problem enters explicitly through the constant C in the energy equation (4): 

C = L c ap\- 2a \v in \ 2a -\ (5) 

where L c sets the reference length scale and a has already been introduced in equation (1). In the follow- 
ing, lengths will be conveniently normalized to the stationary thickness of the cooling region, so that the 
equilibrium position of the shock is x = 1. The explicit expression for C is given in §3.1. 

Relations between quantities ahead and behind the shock follow from the Rankinc-Hugoniot jump 
conditions: 

1 _ 7 - 1 2 

27-1 

Ps ~~ 7TI ~ 7(7 + l)M 2 ' (?) 

where quantities immediately behind the shock are denoted with the subscript s, and 7 = 5/3 will be used 
in what follows. 



3. LINEAR THEORY 



Equilibrium configurations of radiative shock waves may be thermally unstable. The nature of the 
instability may be interpreted as follows. 
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Consider a stationary shock, initially in equilibrium; if, say, the postshock temperature is slightly in- 
creased, a longer cooling path will results and the excess pressure will force the shock to move upward. In the 
frame of the shock, the velocity of the incoming gas will increase even further and the postshock temperature 
will rise according to the square of the preshock velocity. If radiative losses are described by a decreasing 
function of the temperature, the cooling time will increase and the shock will continue to move upward. 



3.1. Perturbative Analysis 

A perturbative study is carried out by properly linearizing equations (2)-(4) around the steady-state 
solutions denoted with p , v and p . The perturbed location of the shock front is written as 

x s = 1 + -/ l , (8) 
o 

where x s = 1 is the shock equilibrium position in absence of perturbation (e = 0), e is the magnitude of the 
perturbation, and S is a complex eigenfrequency. According to the normalization scales introduced in §2, 
time is expressed in units of L c /\vi n \. 

Following the same notations as in Saxton et al. (1998), it is convenient to write 5 = 5^, + iSi, where the 
real part <5r gives the growth/decay term, while Si represents the oscillation frequency. The nature of the 
instability is determined by the sign of <5 R : modes with negative <5r are stable, while modes with positive (5r 
are unstable. 

Perturbed physical variables take the form 

q(t,t) = q (t)(l + \ q (0^ St ) , (9) 

where q G {p, v,p}, go(0 is the corresponding steady-state value and the complex function A q (£) describes 
the effects of the perturbation. Here £ is a spatial coordinate normalized so that £ = 1 at the shock and 
£ = at the wall: 

-J*) (10) 

The fluid equations are linearized in a frame of reference which is co-moving with the shock; in this 
frame the derivatives of a flow variable become 

d d dt; d d dt; d 

di ~* di + lTtirc fa ^ dx~dt;' ( ' 

Therefore, retaining only terms up to first order in e, one has 

|f « (qo^S-tq' Q )ee s \ (12) 

|f - «'o + (l'o\ + 1oX' q |) ee* , (13) 

where a primed quantity denotes a derivative with respect to £. 

The steady-state equations are obtained by collecting the zeroth order terms in the Euler equations; 
conservation of mass and momentum is trivially expressed by 

p Q v = -l, (14) 
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-v +Po = m, (15) 

where the integration constants on the right hand sides may be evaluated from the preshock values; hence 
m= 1 + 1/( 7 A4 2 ). 

The pressure equation provides the explicit dependence on the spatial coordinate; it can be put in closed 
integral form by writing 

where 

/(«)= f {-y?- a [y V {m+ , v)] dy> (17) 

Jo {m + y) a y ' 

and v s = — (l + 3/A / ( 2 )/4 is the fluid velocity immediately behind the shock (eq. [6]). Notice that, according 
to the normalization units introduced in §2, the constant C in equation (5) takes the value 

c = -eFi)- (18) 

Results pertinent to this section are evaluated in the strong shock limit, M — > oo, so m = 1, v s = — 1/4 
and a becomes the only free parameter in the problem. 

The integral in equation (17) can be evaluated analytically for some specific values of the cooling index 
a (CI) but it has to be computed by numerical quadrature in the general case. Notice that a steady-state 
solution is possible only if the integral converges, that is, if a < 3. Equation (16) can be inverted numerically 
to express the postshock steady flow velocity vq as a function of £. The steady-state profiles are shown in 
Figure 1. 

The first-order terms in e provide three coupled complex differential equations for the perturbations; 
using the unperturbed postshock velocity vq as the independent variable, they are 



d\ p d\ v £ X p S d£ 
dvo dv n vf. v dvo 



+ -^ = --2--^^, (19) 



VQ pL ^ po pL = _ K6 *£. + ^ + Xp-2X v -X p , (20) 
dvo dvo dvo vq 



voPo 7 



dv dv 



(2 - a)X p + (a- 1)X P - X v + - 



PoX p 5-^+t, (21) 
dv 



where d£/dvo is given by straightforward differentiation of equation (16) together with equation (17). 

For a given value of a, equations (19) through (21) have to be solved by integrating from the shock 
front, where vq — v s , to the wall, where vo = 0. The eigenmodes of the system are determined by imposing 
appropriate boundary conditions to select the physically relevant solutions. At the shock front the jump 
conditions for a strong shock (M — ► oo) apply (Imamura et al. 1996; Saxton et al. 1998): 

X p = , X v = -3 , X p = 2 . (22) 



At the bottom of the postshock region (£ = 0) the relevant physical solutions must satisfy the stationary 
wall condition, namely, that the flow comes to rest and the velocity must be oscillation-free. This requires 
that both the real and imaginary parts of X v (vo) vanish at v n = 0. The complex frequencies S for which such 
solutions are possible identify the eigenmodes of the system. 
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The method of solution adopted here consists in minimizing the real function of two variables r](8n, Si) — 
|A„(0)|. Here |A„(0)| = (A^ R (0) + A^ j(O)) 1 2 is the magnitude of the velocity perturbation at the bottom of 
the postshock region; the values of the real and imaginary parts, A„ jR (0) and A„ i i(0), are obtained by direct 
numerical integration of equations (19)-(21) for a given pair (<5 R , Si). In practice, since the system is singular 
at the origin, integration proceeds from the shock up to some small value of v , denoted with v e . Setting 
v £ < 10~ 3 did not produce significant variations in the solution. 

A preliminary coarse search with trial values of (5 R and Si shows that, for a given value of a, an 
indefinitely long series of modes exists. Following CI, modes are labeled by increasing oscillation frequency, 
so that n — correspond to the fundamental mode, n = 1 to the first overtone, n = 2 to the second overtone, 
and so on. The approximate position of each mode n, {S^\ S^), is then iteratively improved by repeating 
the search on finer sub-grids (in the complex S plane) centered around the most recent iteration of 
S[ n \ The process stops once the relative error between two consecutive iterations falls below 10~ 6 . For 
practical reasons, the search algorithm has been restricted to the first eight modes for values of a uniformly 
distributed in the range — 2 < a < 2. Results are shown in Figures 2 and 3, while modes for some specific 
values of a are listed in Tables 1 and 2. 

A mode is stable if the real part of the corresponding eigenvalue has negative sign, and unstable other- 
wise. High-frequency modes are characterized by growth rates which decrease faster than low-frequency ones 
for increasing a. Hence, the fundamental mode (n = 0) has the smallest growth/damping rate for a < 1, but 
the smallest damping rate for a > 1. Modes with n > 1 have monotonically decreasing oscillation frequencies 
while, for the fundamental mode, Si reaches a maximum value at a « 1.1 and decreases afterwards. 



3.2. Critical a 

For each mode n, a critical value of the cooling index, ai n \ may be defined, such that 8^ = when 
a = ai n \ Hence, the n-th mode is stable for a > a'"' and unstable when a < etc™' (Fig. 4). The value 
of the critical a is computed by interpolating a with a quartic polynomial passing through the two pairs of 
values across which <5 R changes sign. Thus the fundamental mode becomes stable for a > 0.388, the first 
harmonic for a > 0.782, and so on. Values of a c are listed in Table 2 and shown in Figure 4. Interestingly, 
the sequence of critical a is not monotonic with increasing n. Finally, notice that all (8) modes become 
eventually stable for a > 0.92. 



3.3. Linear Fit 

By inspecting Figure 3, one can easily see that, for a given a, the oscillation frequencies of the different 
modes are approximately equally spaced as n increases. In this respect, they resembles the modal frequencies 
of a pipe open at one end (Toth & Drainc 1993; Saxton et al. 1998; Saxton 1999) and can be described by 
a simple linear fit of the form 

5 (n) = ~ 5 (0) + ^ (23) 

with a small residual, < 0.5%. In equation (23), 5^ is the "fitted" fundamental frequency and A<5i is 
a frequency spacing depending on the cooling index a. A8\ is monotonically decreasing for increasing a. 
Values of 5j (0) and ASi are given in Tables 1 and 2. 
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4. TIME-DEPENDENT NUMERICAL SIMULATIONS 

The results of the previous sections indicate that radiative shocks in real astrophysical settings may be 
linearly unstable and thus far from an equilibrium configuration. This calls for the investigation of the full 
time-dependent problem where nonlinear effects may play a major role in the shock dynamics. 

In what follows, the radiative shock evolution is analyzed through a set of numerical simulations for 
different values of the cooling index a. The early evolutionary phases are of particular interest since they 
can be directly compared to the expected linear behavior, thereby providing an effective tool in validating 
the correctness of the numerical method and choice of boundary conditions. Nonlinear effects, on the other 
hand, describe the long-term dynamics of the shock and play a crucial role in determining whether a linearly 
stable mode may actually be nonlincarly unstable (Saxton 1999). 



4.1. Numerical Method 

The numerical approach followed here relies on the high-resolution shock-capturing methods (HRSC 
henceforth, see LeVeque (1998) and references therein). HRSC methods rely on a finite- volume, conservative 
discretization of the Euler equations, thus being particularly suitable in describing shock dynamics and, in 
general, modeling flow discontinuities. 

The starting point is the system of equations (2) through (4) written in conservative form: 

<9U dF o 

where U = (p, pv, is a vector of conservative quantities, while 



/ 


pv 


\ 






^ 




pv 2 + p 




, s = 







V 


(E+p)v 


) 






-c P 2 - a p a J 



(25) 



are the flux and source term vectors, respectively. The total energy density E is expressed as the sum of 
internal and kinetic terms: 

E 



P pv 
7-1 2 



(26) 



The system of equations (24) is solved numerically using PLUTO, a high- resolution Godunov-type code 
for astrophysical fluid dynamics (Mignone et al. 2005, in preparation). 

With PLUTO, equations (24) are solved by operator splitting, i.e., by treating the advection term dF/dx 
and the source term S in separate steps. This approach is second order accurate in time if the two operators 
have at least the same accuracy and the order into which they are applied is reversed every time step (Strang 
1968). 

During the advection step, a high-resolution, shock capturing Godunov-type formulation is adopted. 
Second-order accuracy in space is based on a conservative, monotonic interpolation of the characteristic 
fields (Colella 1990). Third-order temporal accuracy is achieved by a multi-step Runge-Kutta TVD algorithm 
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(Gottlieb & Shu 1998): 

U] = U? + At n R J -(U n ). 



U ? = i( 3U " + U ) + A ^( Ul )) - (27) 
U™ +1 = 1 (u™ + 2U| + 2Ai"R 7 (U 2 



3 

Here Rj(U) is a conservative, discretized approximation to the flux term on the right hand side of equation 
(24): 

F,.i-F, i 

«-j(U) = - 3+2 3 2 , (28) 

where j labels the computational zone with mesh spacing Axj. The numerical fluxes Fj±i are computed 
using the approximate Riemann solver of Roe (1981). It should be mentioned that, although different 
numerical schemes have also been employed, no significant differences were found in the results presented in 
§4.3. The particular choice of Riemann solver and interpolation algorithm is quite common and represents 
a good trade-off between accuracy and computational time. 

Cooling is treated in a separate source step, where 

*£ = -Ci?-°p°, (29) 

with C given by equation (18), is solved. Notice that only the internal energy is affected during the source 
step, while density and velocity remain constant with the values provided by the most recent advection step. 
Thus, the kinetic contribution to the total energy can be discarded and equation (29) provides an ordinary 
differential equation in the pressure variable only. Integration for a time step At™ can be done analytically: 

p „ +1 = < | [p 1 - a -At"C( 7 -l)/> 2 - a (l-a)]^ if a^l, ^ 
pexp(-C( 7 - l)AiV~") if a = l, 

where p and p are the pressure and density at the beginning of the source step and the suffix j has been 
omitted for clarity of exposition. 

Radiative losses are identically zero for T < T c , where T = p/p is a dimensionless temperature and T c 
is the cutoff temperature, equal to the temperature of the incoming gas, i.e., T c = 1/(^M 2 ). 



4.2. Initial and Boundary conditions 

The computational domain is the region xo < x < x\, where Xo and x\ define the locations of the lower 
and upper boundaries, respectively. As in the previous sections, lengths are expressed in units of the cooling 
thickness in steady-state; hence the shock is initially located at x = 1. 

Upstream, for 1 < x < x\, density, velocity and pressure are uniform and equal to the prcshock values, 
i.e. p = 1, v = — 1 and p = l/(jM 2 ). At x = x\ a constant supersonic inflow defines the upstream boundary 
condition. Here x\ = 5 and M. = 40 will be used for all simulations. Downstream, for xo < x < 1, flow 
quantities are initialized to the steady-state solution given by equations (14) through (16). The location of 
the lower boundary xq depends on the particular boundary condition adopted. 
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The implementation of a suitable boundary condition at the lower boundary poses a serious nontrivial 
problem for a number of different reasons. Close to the wall, density and velocity experience rapidly varying 
sharp gradients, thus demanding increasingly high resolution in order to resolve flow patterns. This, in turns, 
has the inevitable consequence of restricting the time step size in an explicit code. Although a number of 
different strategies have been proposed, there is no general agreement about what would be a "consistent" 
boundary condition. A comparison between four different approaches is considered in this paper and results 
are discussed in detail in the next sections. 

A first approach consists in using a reflecting boundary condition (RFbc), commonly adopted (Plewa 
1995, SB) to simulate the presence of a rigid wall or to enforce axial or planar symmetry. It imposes 
symmetric profiles (with respect to the "wall" position xo — 0) on density and pressure, while the velocity 
is antisymmetric, i.e., v(—x) = —v(x). Hence at the lower boundary, the velocity is zero at all times while 
density and pressure have zero gradient. 

A second approach adopts a "fixed" , time- independent boundary conditions (FXbc) where flow variables 
are kept constant at the steady-state values at x — x a . In this case, x a = 10~ 3 can be safely used, since is put 
sufficiently close to the wall but such that the temperature at that point is still above the cutoff temperature. 
The cutoff temperature, therefore, has only the effect of preventing cooling in the upstream region. 

I also introduce a new, third approach based on the characteristic boundary method (Thompson 1987, 
1990), the details of which are given in Appendix A. Following the same approach used in §3.1, the velocity 
will be held constant at the steady-state value, whereas pressure and density are allowed to evolve with time. 
The "constant-velocity" boundary condition is hereafter referred to as CVbc. 

Finally, a fourth recipe (SB, SDB, LCSb) is to further extend the domain in the downstream direction 
by placing a cold dense layer for xi ow < x < xq. Here xi ow = —2 is the new location of the lower boundary, 
while x — 10~ 4 still lies inside the postshock region. Flow quantities have constant profiles, continuous at 
x = x and the temperature of the cold layer becomes approximately 2% of the temperature immediately 
behind the shock. For this boundary condition, the cutoff temperature was set equal to the temperature of 
the cold dense layer. At the back of the layer, x — x\ ov/ , an outflow boundary condition (OBbc) holds on 
density, velocity and pressure. 

The onset of instability is triggered by the discretization error of the numerical scheme and no external, 
"ad hoc" perturbation is introduced, unless otherwise stated. 

A uniform grid is used in the region xo < x < 1.4, whereas a second, geometrically stretched grid covers 
the rest of the upstream region, 1.4 < x < 5. The extent of the uniform region has been chosen to ensure 
that the largest-amplitude oscillations would adequately be resolved. 

Issues concerning grid resolution effects must not be underestimated. In fact, as outlined by Sutherland 
ct al. (2003, SBD hereafter), sharp density gradients can be described with relatively limited accuracy 
because of numerical diffusion effects that cause high-density regions to "leak" mass into neighboring low- 
density zones. Since the cooling process is proportional to the square of density, radiative losses will generally 
be overestimated, causing abnormal, excessive cooling. Although this issue is intrinsic to any grid of finite size 
and cannot be completely removed, higher resolutions can considerably mitigate the problem. Furthermore, 
small-amplitude oscillations of the shock front can be adequately captured on finer grids. 

For this reason, a grid of 2240 computational zones covers the extent of the uniform region for all 
numerical simulations presented in this paper. With this resolution, the postshock flow is initially resolved 
on 1600 points. The number of points for the stretched grid has virtually no influence on the solution and 
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is held fixed at 200 in all cases. 



4.3. Results 

The one-dimensional numerical simulations are carried out for five different values of a, selected accord- 
ing to the linear analysis results presented in §3. 

I first consider the case where a = 0, since all modes have positive (linear) growth rates. Next, a = 0.5 
and a — 0.7 are examined, since all the overtones save the fundamental have positive linear growth rates 
and are, therefore, unstable. The value a = 0.8 lies just above the instability limit of the first and second 
harmonics; therefore, as one can see from Table 2, only overtones with n > 3 are expected to be unstable. 
Finally, I consider the case a = 1, for which all the first eight modes have been shown to have negative 
growth rates and thus are stable. 

4.3.1. a = 0.0 

The time-evolution diagrams for the four boundary conditions presented in §4.2 are shown in the four 
panels of Figure 5 and the corresponding power spectra of the shock oscillations are shown in Figure 6. 

Significant departure from the steady-state solution occurs most rapidly when the RFbc is employed. 
In this case, the amplitude of the oscillations rapidly increases with time and the system enters a nonlinear 
saturated regime around t ~ 20 after a short-lived linear phase. The use of the RFbc yields oscillation 
frequencies which are found to be shifted with respect to the values obtained from linear analysis. This 
particular choice of boundary condition, in fact, forces the velocity to have a node at the location of the 
lower boundary, while density and pressure have an extremum. Since this condition is far from the equilibrium 
configuration (cqns. [14]-[16]), strong nonlinear perturbations originate in the postshock flow and steepen 
into secondary shocks at a high rate. This also contributes to the higher amplitude oscillations observed in 
this case. Not surprisingly, the power spectrum of the shock position, Figure 6, exhibits frequencies that 
are offset from the ones predicted by linear analysis. Notice that, since the mass flux through the lower 
boundary is zero and cooling is not effective for T < T c , a cold layer of gas at T = T c accumulates at the 
bottom of the cooling region (Fig. 5). Inside this layer, density is approximately constant and equal to -fM 2 , 
whereas waves propagate at the local sound speed, c s ~ 1/A4. 

In contrast, the CVbc and FXbc yield similar results and the system preserves profiles close to the 
initial steady-state values. The early phase of evolution (t < 60 ~ 70) is characterized by a linear growth of 
the perturbation, while, for t > 80, the amplitude of the oscillations begins to saturate and the instability 
becomes nonlinear. During this phase the largest oscillation amplitudes reach ~ 25% of the initial equilibrium 
position. The OBbc shows reduced amplitudes with respect to the previous cases. In addition, the linear 
phase is longer than the CVbc or FXbc, and the transition to the nonlinear regime occurs only for t> 110. 

The power spectra for the early phase of the evolution (0 < t < 41) are plotted in Figure 6. Both the 
CVbc and FXbc yield eigenfrequcncies that can be definitely identified with the theoretical values, with a 
bigger uncertainty in the fundamental mode, see Table 3. Results obtained with the OBbc are similar and 
the fundamental mode differs from the analytical expectation by less than 4%. 

One should bear in mind, however, that the linear growth rate of the fundamental mode is a factor of 
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6 smaller than those of higher harmonics (1 < n < 7, see Table 2), and therefore modes with n > 1 tend 
to saturate faster. For this reason, a representative sample of the linear phase must have a limited length 
in the time domain and, consequently, the power spectrum inevitably suffers from poor resolution at lower 
frequencies. 

The situation is different when a longer portion (92 < t < 201) of the shock position during the saturated 
regime is analyzed, see Figure 7. In this case, the predominant mode of instability is the first overtone, 
whereas the fundamental mode and the second harmonic contribute by less than 10% to the oscillatory 
cycle. Notice also that when the CVbc and FXbc are used, the prevalent frequency of oscillation differs from 
the linear prediction by ~ 10%, but it coincides with the first harmonic when the OBbc is adopted. 

In all cases, a main sequence of harmonics with increasing frequencies fl^\ Q( u \ Q,( m \ etc., may 
be identified by inspecting Figure 7. These overtones have monotonically decreasing power and do not 
necessarily coincide with the linear modes, but result from nonlinear interactions. Nonlinearity enters through 
mode- mode coupling, as it is suggested by considering the frequency spacing between them (Table 4) . For the 
CVbc and OBbc, in fact, the spacing appears to be cither a multiple of the fundamental mode or equal to the 
first overtone. In the OBbc case, for instance, one has that 0( n )-fiW w O^ 111 ' w 0< IV ) -0( m ' « fi« , 

i.e., the mode spacing is a multiple of the first overtone. A similar result holds in the CVbc case where, from 
Table 4, it can be verified that Q( n ) - fi« rj 3ft( \ Q (u ^ - ft< n ) « ft< IV ) - 0( m ) « Q«, and so on. 

Intermediate, secondary peaks associated with small-power modes appears between the main sequence 
overtones. Some of these modes have been identified and labeled in Figure 7withft( Ia \ ft( Ib \ tt( IIa \ ft( IIb ), 
etc.. These secondary overtones may result from mode- mode coupling between the main sequence modes 
and the fundamental. This coupling is most evident for the CVbc, where one finds that £l( 0a ) — w 
ft(ia) _ ~ ft(ii) _ fi(ib) „ q(0) ) and S i m ii ar i y f or higher harmonics (sec Tab. 4). 

4.3.2. a = 0.5, 0.7 

Based on the previous results and considerations, the CVbc, FXbc and OBbc yield results that more 
accurately reproduce the predicted linear behavior during the system's early phase of evolution. Moreover, 
results obtained with the FXbc and CVbc exhibits strong similarities, and thus only the CVbc and OBbc 
will be considered in what follows. 

The a = 0.5 value is of particular astrophysical relevance, since it describes optically thin bremsstrahlung, 
which is the main source of radiative losses at temperatures of the order of 10 8 -j- 10 9 K, typical in accretion 
shocks in magnetic cataclysmic variables. 

When the CVbc is adopted the system exhibits a linear phase for t < 150, gradually followed by the 
transition to the nonlinear regime. When compared to the a = case, the oscillation amplitudes in the 
saturated regime are reduced by a factor of approximately 50%. The situation is quite different, however, 
when the OBbc is considered: Figure (8) shows that the solution remains close to the initial steady-state 
values and unstable oscillations grow at a smaller rate. 

A similar behavior has been reported by SB95 (who also adopted an "open boundary" condition) for 
small values of the inflow Mach number (i.e., M. — 5). In their simulations, however, the amplitude of 
the oscillations was found to increase for higher Mach numbers, a behavior not observed in the present 
work. The present conclusion is supported by several supplementary tests in which both the inflow Mach 
number and the density of the cold layer were changed, but a fully nonlinear growth of the instability was 
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still never observed. In all the numerical tests, in fact, the cold gas layer always acts as an absorber to 
incoming perturbations, consequently reducing the amplitude of the reflected waves. Even in the presence 
of an external "ad hoc" perturbation (similar to the one introduced in SB95) it was found that the use of a 
cold dense "layer" inhibits the growth of instability when a > 0.45. 

The behavior of the system during the early phases is reflected in the power spectra shown in Figure 9, 
where a positive identification of the oscillation eigcnfrequencies with the linear ones is clear. The relative 
errors of the identifiable peaks are less than 4% for all modes, see Table 3. Notice that the fundamental mode 
(expected to be stable from the linear analysis) is also visible in the spectrum, since the initial numerical 
perturbation excites all modes regardless of their stability. 

The power spectra taken during the later phases, Figure 10, reveal that, for the OBbc, the (mildly) 
unstable behavior is mostly sustained by the first three overtones, whereas the first harmonic is the only 
dominant mode for the CVbc. In both cases, little contribution is given by the fundamental mode. Nonlinear 
effects, however, suggest that the fundamental mode may still be important through mode-mode coupling. 
Similarly to the a = case, in fact, a main sequence of modes can again be identified, see Figure 10. For the 
CVbc, the frequency spacing between these modes is either a multiple of the fundamental or equal to the first 
overtone, e.g., ft( n ) - fiW w 3fi(°) and fi< IV ) - O^ 111 ' w fiW. Secondary, small -power overtones are mainly 
visible for the OBbc case. Again, strong evidence for inter-mode coupling is supported by the fact that these 
secondary overtones may be decomposed into main sequence modes. In fact, if one consider the frequencies 



listed in Table 4, it can be seen that fl^ « Q« - fi( n ) - 0< Ia ) w fi( IIIa ) - fi( m ) w and so 



For a — 0.7, an additional external perturbation has been introduced to catalyze the onset of instability. 
The perturbation is initially given in the velocity profile as 



with e = 0.05. Density and pressure are obtained according to equations (14) and (15). 

The different behaviors of the OBbc and CVbc are illustrated in Figure 8. Results obtained with the 
OBbc show that the initial perturbation is damped roughly on a timescale t ~ 150. As for the a — 0.5 case, 
the cold dense layer behind the postshock region tends to quench large-amplitude perturbations. On the 
other hand, when the CVbc is adopted, the initial perturbation does not fade away and the instability grows 
at a small rate. The amplitude of the oscillations relative to initial shock position is now further reduced to 
< 5% of the initial shock position. Table 3 shows that the eigenfrequencies of the oscillations differ by less 
than 8% from the theoretical results. 

The power spectra for the early linear phase, Figure 9, are similar to the previous cases, although only 
modes with < n < 4 (for the CVbc) and 1 < n < 5 (for the OBbc) contribute to the oscillations. 

During the nonlinear phase (CVbc only), the first harmonic gives the largest contribution, while the 
third overtone account for roughly 10%, Figure 10. Although little power is present in the fundamental mode, 
the frequency spacing between main sequence harmonics seems to indicate that mode-mode coupling may 
account, one more time, for the secondary, small-power intermediate peaks (fl^ , ^ Ih \ etc.) shown in Figure 



10. Some the O's are, in fact, closely related: ^ w fiW-fiW, n (Ic 1 « Q( m ) -Q( Ia ), ft( IIIa ) - ft( m ) w , 

so on. 



on. 




(31) 
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4-3.3. a = 0.8,1 

The simulations for the last two cases, a — 0.8 and a = 1, are carried out using the CVbc only, since 
no growth of instability was observed using the OBbc. In both cases, the perturbation given by equation 
(31) was imposed at t = 0. Results are shown in Figures 11 and 12. 

When a = 0.8, the early phases of the evolution reflect the expected linear growth, as one can see 
from Figure 12. The power spectrum for to this phase shows modes of oscillations that clearly match the 
theoretical ones. Most of the power is contributed by the first harmonic, followed by the fundamental and 
then the remaining overtones. 

The complete transition to the nonlinear phase occurs for t > 200, where smaller-amplitude, higher- 
frequency oscillations take over. A power spectrum of the late evolution reveals the effects of this transition, 
Figure 12. Most of the power is concentrated in the third harmonic, with only ~ 10% going into the first 
overtone and less than ~ 1% into the second harmonic. The fundamental mode is absent. 

The fact that high-frequency oscillations arc dominated by the third harmonic is quite a remarkable 
result, since this overtone is the lowest unstable mode only in the narrow range 0.795 < a < 0.85, while 
modes with n < 3 are stable as can be seen from Table 2. This strongly suggests that this particular choice 
of boundary condition is particularly consistent with linear results. 

Notice that both the first and second harmonic should be linearly stable, since = 0.7815 < 0.8 and 
af = 0.795 < 0.8. Their presence in the power spectrum, however, indicate that a weak nonlinearity may 
probably be present. Besides, nonlinear interactions are likely to be responsible for the frequency coupling 
between the third and fourth main sequence overtones, since f)( IV ) w 2Q( m ). 

Finally, when a = 1, the initial perturbation is damped and the system returns to the original equilibrium 
solution for t > 70. The power spectrum of the early evolutionary phase (Fig. 12) shows that the shock 
oscillations are decomposed into frequency modes that are well approximated by linear results. 



5. DISCUSSION 

A study of planar radiative shocks with a power-law cooling function A ~ p 2 T a has been conducted. 
Both linear and nonlinear time-dependent calculations have been presented. 

A linear stability analysis has been carried out for several values of the cooling index a in the range 
[—2, 2). For a given value of a, multiple discrete modes of oscillation exist and the real and imaginary parts 
of the first eight cigcnfrcqucncies have been derived. The overstable modes are labeled in order of increasing 
oscillation frequency so that n = corresponds to the fundamental mode, n = 1 to the first overtone, 
n = 2 to the second overtone, and so forth. The stability criterion of a particular mode is expressed by 
the condition a > ai n \ where a[™' is the critical value of the cooling index for the n-th mode. For the 
fundamental mode, for example, = 0.388, whereas for the first and second harmonic = 0.782 and 

(2) 

a c = 0.795, respectively. A general trend towards stability exists for increasing a, so eventually all modes 
are stabilized for a > 0.92. This study confirms previous results (CI), for which only the first two or three 
modes have been reported for some values of the cooling index. It has been shown that oscillation frequencies 
are linearly proportional to the mode number n, a behavior similar to the quantized modes in a pipe. 

The perturbative study has been complemented by several numerical simulations using an Eulerian, high- 
resolution shock-capturing scheme. The shock evolution has been followed through the linear and nonlinear 
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phases for different values of a and boundary conditions. Among the four boundary conditions under 
consideration (§4.2), a new time-dependent boundary treatment has been introduced. The new approach is 
based on the characteristic boundary method for the Euler equations and is particularly consistent with the 
basic assumptions used in the analytical work, where the velocity perturbation has a node at the wall. 

For the most unstable case considered here (a — 0), all boundary conditions yield similar results, 
although the reflective wall boundary condition is not particularly suitable in modeling small departures 
from the stationary solution. The remaining three strategies provide modes of oscillations that, in the limit 
of small perturbations, are close (within 5% accuracy with the exception of the fundamental mode) to the 
analytical values. 

The cases a = 0.5,0.7,0.8 and 1 have also been considered. For a = 0.5,0.7, the numerical simulations 
show that the choice of the lower boundary condition has a more severe impact on the growth of unstable 
modes during the saturated phase. For example, the additional cold dense layer used in the open boundary 
condition (§4.2) tends to inhibit large-amplitude oscillations and, for a > 0.7, totally prevents the growth 
of instability. In contrast, the new boundary approach yields results that closely reflect the analytical 
predictions; unstable behavior was observed for a = 0.5, 0.7, and 0.8, although the saturated amplitude of 
oscillations considerably decreases with increasing a. For a = 0.8, the largest oscillations during the saturated 
phase are reduced to ~ 0.5% of the initial shock position and are thus barely visible at the resolution adopted 
(1600 zones for the postshock flow). For a = 1, the shock is stable and initial perturbations are damped 
on a characteristic timescale roughly proportional to the e-folding time of the first overtone. The modes 
of oscillations found in the numerical simulations (during the early phase of evolution) can be positively 
identified with the ones derived by linear analysis. The relative error is usually small, < 8%. Notice that 
this error also accounts for the discrete frequency spacing introduced by the fast-Fourier transform of the 
shock position. 

The transition from the linear to the nonlinear regime has also been investigated. Power spectra of the 
shock position during the late evolutionary phases reveal that the first overtone is the dominant mode of 
oscillation when a < 0.7, but that the third harmonic contributes to most of the power at a = 0.8. The 
contribution of the fundamental is only 10% for the most unstable case (a = 0), and decreases for increasing 
a. 

The new result of this work shows that a main sequence of overtones characterizes the saturated, 
nonlinear oscillatory phases. Additional, secondary modes may also be present, depending on the particular 
choice of boundary condition. These modes of oscillation do not always match those predicted by linear 
analysis but result from complex nonlinear interactions. For the first time, evidence has been provided in 
favour of mode-mode coupling, particularly between the first harmonic and the fundamental mode. The 
result extends also to those cases in which some of the modes are linearly stable (a — 0.5,0.7,0.8), thus 
supporting the possibility that linearly stable modes may actually become nonlinearly unstable. 

In summary, a general trend towards stability is found for a > 0.8, while an unstable behavior is expected 
for a < 0.4, regardless of the choice of the lower boundary condition. On the other hand, numerical models 
of radiative shocks are more sensitive to the treatment of the lower boundary condition when 0.4 < a < 
0.8, a range particularly relevant when optically thin bremsstrahlung is the dominant cooling mechanism. 
It should be pointed out that the use of a cold layer of finite thickness may be more self-consistent in 
realistic astrophysical applications. The existence of the layer is, in fact, automatically induced by a cutoff 
temperature in the cooling function and avoids the complication of specifying a boundary conditions at the 
interface between the postshock flow and the layer (provided the sharp density gradients present in this 
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region are adequately captured). 

In spite of the oversimplifying assumptions adopted in this study these results show a number of 
interesting consequences for a variety of astrophysical settings. 

Radiative shocks with velocities v s > 130 km s -1 are not uncommon in jets from young stellar objects, 
supernova remnants in the radiative phase, magnetospheric accretion in T-Tauri stars, and colliding stellar 
winds in relatively close binary systems. For these systems, the shocked interstellar gas reaches temperatures 
in the range 10 5 — 10 7 K and cools mainly by line emission, for which a < —0.5. Under these conditions, 
radiative shocks are likely to show unstable behavior in all modes and phenomcnological interpretations 
based on steady-state models become of questionable validity (Innes et al. 1987a, b). Although inclusion 
of transverse magnetic fields extends the range of stability (Smith 1989; Toth & Draine 1993), the global 
thermal instability of radiative shock waves may still be important in interpreting a number of distinct 
observational features, such as emission-line ratios observed in interstellar radiative shocks (Hartigan et al. 
1994), mixing between hot and cold material in colliding winds (Stevens et al. 1992; Antokhin et al. 2004), 
the filamentary structures observed in supernova remnants (Blondin et al. 1998; Walder & Folini 1998), and 
so forth. 

Less conclusive assertions can be made for standing shocks in the accretion columns of Polar and 
Intermediate Polar systems. At temperatures of the order of 10 8 — 10 9 K the X-ray emission is primarily 
determined by optically thin bremsstrahlung, although cyclotron and Compton cooling may not be neglected 
(Saxton et al. 1998). However, in the simple case where radiative losses are due to bremsstrahlung cooling 
only, a « 0.5, the dynamics of the shock may be influenced by the interaction with the upper photosphcric 
layers of the white dwarf (Cropper 1990). Hence realistic models of accretion columns may require a more 
complex treatment of the lower boundary. For this reason, inclusion of additional physical processes such 
as magnetic fields, multi-dimensional effects, thermal conduction, etc., might be crucial for drawing firm 
conclusions about the stability of radiative shocks in AM-Her-type systems. Some of these issues will be 
considered in future extensions of this work. 

The author would like to thank T. Plewa, B. Rosner, A. Konigl, D. Q. Lamb, and T.J. Linde for their 
constructive support during my Ph.D. research program at the University of Chicago. 

A. CHARACTERISTIC BOUNDARY CONDITIONS 

The hyperbolic nature of the Euler equations requires boundary conditions to be specified according 
to the way information propagates in and out of the boundary. The novel approach introduced in §4.2 is 
based on the characteristic boundary method (Thompson 1987, 1990), where "physical" and "numerical" 
boundary conditions specify how zone boundary values are integrated in time along with the interior values. 
Although the subject of boundary conditions is a vast one and falls outside the scope of this paper, details 
of implementation are given hereafter. 

A "physical" boundary condition describes information that enters the domain and thus has to be 
imposed for each characteristic wave that propagates from the boundary towards the inside. Information 
directed outside the boundary is entirely determined by the solution inside the domain and thus does not 
require a boundary condition. The numerical scheme, however, still depends on the knowledge of all flow 
variables at boundary zones, and hence additional "numerical" boundary conditions must be prescribed in 
a consistent way. 
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In the present context, the boundary equations are more conveniently formulated using the quasi-linear 



form 
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where r = 1/ p. The system of equations (Al) holds at the boundary and must evolve in time together with 
the interior values. The characteristic speeds of the system (Al) are given by A ± = v ± c s , A = v, where 
c s = ("ip/ ' p) 1 ! 2 is the speed of sound. 

In the postshock region the flow is initially subsonic everywhere, since — c s < v < for < x < 1. In 
the limit of small perturbations around the steady-state values, it is reasonable to assume that a condition 
for subsonic outflow will continue to hold at subsequent times. Hence, the characteristic associated with A + 
has positive sign (i.e., it carries information inside the domain), whereas A and A~ are directed outward. 
This means that only one "physical" boundary condition can be freely specified (e.g., a constant pressure 
or velocity) and the remaining two must be compatible with the interior discretization scheme. Choosing a 
constant outflow velocity, for example, is consistent with the linear perturbative analysis, where the velocity 
perturbation is forced to have a node at the origin. 

Integration of the boundary equations (Al) proceeds by splitting the time-dependent solution into a 
contribution coming from the steady-state value and a time-dependent "deviation" : 

q(x,t)= qi (x,t)+qo(x) (A2) 

where q £ {r, u,p}. In the "constant- velocity" boundary condition, the velocity perturbation v\ — at all 
times and the boundary equations prescribe how pressure and density should evolve with time: 

— + v — - t— - t — (A3) 
dt dx dx dx 

dpi dpi dv x dv Q 

-dT + v ^ + "to = ^ - 1)(A - Ao) (A4) 

where the spatial derivatives are computed using one-sided approximations. Notice that equations (A3) and 
(A4) are not a linearization around a stationary solution but are, in principle, valid for arbitrary departures. 
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Fig. 1. — Steady-state profiles for density, temperature and velocity when a = 0. The "wall" is located at x = and supersonic 
gas flows from the right to the left. Flow variables are normalized to their inflow values, and the abscissa is expressed in units 
of shock height. 




Fig. 2. — Growth/damping rates for the first 8 modes as function of a. The solid line represents the fundamental mode n = 0, 
whereas the different symbols (described by the legend in the upper-right portion of the plot) correspond to the seven overtones 
1 < n < 7. Eigcnmodes with 5r < are stable, whereas modes with <5r > are unstable. 
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Fig. 3. — Oscillation frequencies for the first 8 modes as function of a. The symbols have the same meaning as in Figure 2. 
Modes with 1 < n < 7 have oscillation frequencies that are monotonically decreasing functions of a. For the fundamental mode 
(n = 0), however, <5[ increases to reach a maximum at approximately a = 1.1 and decreases afterwards. 
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Fig. 4. — Critical value of the cooling index as function of the mode number n. For a given mode n, values of a > a c have 
negative growth rates and thus are stable. 
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Fig. 5. — Space-time diagrams for a = and the four different boundary conditions described in the text. From top 
to bottom: the reflective boundary condition (RFbc), the fixed boundary condition (FXbc), the constant velocity boundary 
condition (CVbc) and the "open boundary" condition (OBbc). The spatial coordinate is represented on the vertical axis, 
whereas the time evolution of the system is described by the horizontal axis. Time is expressed in units of L c /\vi n \, where L c 
is the length of the cooling region in steady-state, see §3.1, and v ln is the fluid velocity ahead of the shock. In each panel, 
the gas flows supersonically from the top to the bottom. The gray-scale shows the density logarithm: lighter (darker) shades 
corresponds to lower (higher) density regions. 
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a = 0, RFbc , 0.0 < t < 40.9 
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Fig. 6. — Power spectra of the shock position for a = 0.0 derived from the numerical simulations for the four boundary 
conditions described in the text. On the top: reflective (left) and constant velocity boundary conditions (right); on the bottom: 
the fixed (left) and open boundary condition (right). The vertical axis represents power on a logarithmic scale, normalized to 
its maximum value. The horizontal axis shows frequency on a linear scale. The power spectra are obtained by computing the 
Fourier transform of the function x s h(t) — x s h(0), where x s h(t) is the shock position at time t. The transform is taken over 
a sample of length < t < 40.9 The dashed vertical lines in correspondence of the S„, with < n < 7, mark the oscillation 
cigenfrcqucncics derived from linear analysis, see Table 1. 
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Fig. 7. — Same as Figure 6, but when a longer portion of the shock position at later evolutionary times (92.1 < t < 201) 
is Fourier-transformed. For the CVbc and OBbc, a main sequence of overtones (I, II, III and so forth) with frequencies QW, 
, may be identified from the plots. These modes are almost equally spaced in frequency and may result from 
nonlinear mode-mode coupling between the fundamental mode and the first overtone. Secondary harmonics with oscillation 
frequencies ft< Ia ), n< IIa ), ft( IIb ), etc., appear between the main sequence modes. The explicit values of the H's are given in 
Table 4. 




Time 



Fig. 8. — Evolutionary space-time diagrams for a = 0.5, 0.7 with the OBbc and CVbc. Oscillation amplitudes arc considerably 
reduced when a increases and also when the OBbc is adopted. In the worst case (bottom panel), the initial perturbation is 
damped and the system returns to a stationary, stable configuration. 
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a = 0.5, CVbc , 0.0 < t < 42.4 
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Fig. 9. — Power spectra for a = 0.5 (top) and a = 0.7 (bottom) during the early phases of evolution, < t < 42.4 and 
< t < 43.3, respectively. Results obtained with the constant velocity and open boundary condition are shown on the left and 
on the right, respectively. The dashed vertical lines correspond to the frequencies derived from linear analysis. The vertical 
and horizontal axis have the same meaning as in Fig. 6. 
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a = 0.5, CVbc, 236.8 < t < 350.0 a = 0.5, OBbc, 236.8 < t < 350.0 




Fig. 10. — Same as Fig. 9, but for a longer portion of the oscillatory cycle during the late evolutionary phases, 236.8 < t < 350 
(for a = 0.5) and 234.8 < t < 350 (for a = 0.7). During this time window, the oscillation amplitudes have fully saturated and 
a main sequence of modes, similar to those described in Fig. 7, may be distinguished in the power spectra. The corresponding 
oscillation frequencies are labeled by fiW, fi( n ), f2( m ', and so on. A number of secondary harmonics (f2( Ia ), f2' IIa ', Q' IIb ', 
etc.) is present as well. Notice that the fundamental mode Q(°\ although linearly stable, is still present (with little power) 
in the spectra. Values of the different Q's are listed in Table 4. When a = 0.7 and the OBbc is used, the system returns to 
steady-state and a flat spectrum is obtained. 
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Fig. 11. — Space-time evolutionary diagram for a = 0.8 (on the top) and a = 1 (on the bottom) is shown. Only the CVbc 
has been employed. In order to make the oscillations more visible, the spatial scale in the plot shows a small area around the 
shock position. 
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Fig. 12. — Power spectra for a = 0.8 (top panel) and a = 1 (bottom panel). The earlier (0 < t < 43.9 for a = 0.8 and 
< t < 45.4 for a = 1) and later (233.3 < t < 350 for a = 0.8 and 229.3 < t < 350 for a = 1) simulation phases arc shown 
to the left and to the right, respectively. Only the a = 0.8 case evolves into a (weakly) nonlinear phase, characterized by very 
small-amplitude oscillations. Notice how the first and second harmonics, although linearly stable, are still present during the 
later simulation phases (QW and Q( n ) on the top, right panel). The fundamental mode is practically absent. Explicit values 
of the nonlinear overtones are given in Table 4. When a = 1, the initial perturbation produces modes that can be clearly 
identified with the linear ones (bottom, left panel). These modes are absent during the later phases (bottom, right panel), thus 
confirming that the shock is stable. 



-32 - 



Tabic 1. Real and imaginary parts of the complex eigenfrcquencics S — + iSi are given for the first 
eight modes, n — 0..7, and for negative values of a. The lower portion of the Table shows the coefficients 
derived from the linear fit = 5^ + nASi, where 5^ is the "fitted" fundamental mode and A<5i is the 

frequency spacing. 



a = -2 a = -3/2 a = -1 a = -1/2 



Mode <5r <5i <5r 8j 5r <5i <5r <5i 



n = 





0.1671 


0.2175 


0.1353 


0.2416 


0.1031 


0.2616 


0.0693 


0.2787 


n = 


1 


0.3443 


0.9581 


0.2925 


0.9566 


0.2393 


0.9510 


0.1827 


0.9399 


n = 


2 


0.3905 


1.7252 


0.3360 


1.7052 


0.2786 


1.6778 


0.2161 


1.6398 


n = 


3 


0.4258 


2.4622 


0.3707 


2.4277 


0.3121 


2.3820 


0.2476 


2.3204 


n = 


4 


0.4538 


3.2200 


0.3957 


3.1704 


0.3334 


3.1059 


0.2642 


3.0197 


n = 


5 


0.4684 


3.9594 


0.4110 


3.8945 


0.3495 


3.8112 


0.2812 


3.7012 


n = 


6 


0.4918 


4.7110 


0.4319 


4.6330 


0.3669 


4.5325 


0.2938 


4.3996 


n = 


7 


0.4982 


5.4548 


0.4389 


5.3604 


0.3749 


5.2400 


0.3039 


5.0819 



5(0) 



A<5i 



3°> 



A5i 



5^ 



ASi 



<5<°> 



A<5i 



0.2183 0.7486 



0.2352 0.7324 



0.2502 0.7129 



0.2642 0.6882 



Table 2. Real and imaginary parts of the complex eigenfrequencies 6 = 5n + i5\ are given for the first 
eight modes, n = 0..7, and for nonncgative values of a. The rightmost column gives the critical value of a 
for a given mode n, such that for a > the n-th mode is stable. The lower portion of the Table lists the 

coefficients derived from the linear fit. 



Mode 


a 


= 


a = 


1/2 


a = 


1 


a = 


3/2 




Sr 


Si 


Sr 


Si 


Sr 


<5i 


Sr 


Si 


n = 





0.0323 


0.2934 


-0.0101 


0.3052 


-0.0622 


0.3121 


-0.1346 


0.3054 


0.3881 


n = 


1 


0.1201 


0.9210 


0.0476 


0.8887 


-0.0420 


0.8307 


-0.1668 


0.7075 


0.7815 


n = 


2 


0.1450 


1.5857 


0.0602 


1.5043 


-0.0485 


1.3698 


-0.2020 


1.1022 


0.7949 


n = 


3 


0.1739 


2.2347 


0.0851 


2.1087 


-0.0310 


1.9068 


-0.2141 


1.5186 


0.8822 


n = 


4 


0.1841 


2.9003 


0.0865 


2.7242 


-0.0396 


2.4386 


-0.2128 


1.9042 


0.8578 


n = 


5 


0.2021 


3.5505 


0.1052 


3.3322 


-0.0280 


2.9850 


-0.2416 


2.3021 


0.9111 


n = 


6 


0.2081 


4.2158 


0.1030 


3.9454 


-0.0307 


3.5106 


-0.2390 


2.7186 


0.8959 


n = 


7 


0.2214 


4.8666 


0.1188 


4.5561 


-0.0286 


4.0605 


-0.2375 


3.0974 


0.9196 



A<5, 



5<°> 



A<5r 



A<5r 



A<5, 



0.2774 0.6553 



0.2898 0.6088 



0.3011 0.5359 



0.3076 0.3998 
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Table 3. Relative errors of the oscillation frequencies found from the numerical simulations during the 
early linear phases. The errors are computed as \uj < f l ' > / 5 (n ' ) — 1|, where uj^ corresponds to the closest 
frequency peak referred to the theoretical value. Notice that the finite length of the time window At over 





which the Fourier transform 


is taken 


introduces an uncertainty ~ 1 /At. 




Mode 


a 


= 


a = 


1/2 


a = 


0.7 


a = 0.8 


a = 1 


CVbc 


OBbc 


CVbc 


OBbc 


CVbc 


OBbc 


CVbc 


CVbc 


n = 


0.482 


0.036 


0.039 




0.071 




0.087 


0.122 


n = 1 


0.010 


0.010 


0.008 


0.008 


0.010 


0.010 


0.010 


0.010 


n = 2 


0.042 


0.042 


0.025 


0.025 


0.017 


0.017 


0.012 


0.0003 


n = 3 


0.048 


0.020 


0.026 


0.026 


0.016 


0.016 


0.009 


0.006 


n = A 


0.005 


0.005 


0.031 


0.031 


0.017 


0.017 


0.009 


0.011 


n = 5 


0.016 


0.027 


0.012 


0.032 


0.018 


0.018 


0.010 


0.010 


n = 6 


0.009 


0.009 


0.004 


0.004 


0.018 


0.018 


0.008 




n = 7 


0.032 


0.001 


0.035 


0.035 


0.046 


0.013 


0.010 


0.012 



Table 4. Nonlinear frequencies relative to the later evolutionary phases for a — 0,0.5,0.7 and 0.8. 
Frequencies labeled with f^ 1 ), Q^ u \ fi( in ), and so on, identify the main sequence overtones, whereas the 
intermediate secondary modes are enumerated by appending a letter to the main sequence mode number 
that precedes it (i.e. ft< Ia \ ft( IIa \ fl^ Uh \ etc.) The error introduced by the Fourier transform is ~ 1/ At, 
where At is the time window over which the transform is taken. 



Mode 


a 


= 


a = 


1/2 


a = 0.7 


a = 0.8 


CVbc 


OBbc 


CVbc 


OBbc 


CVbc 


CVbc 


n(°) 


0.287 


0.287 


0.278 


0.278 


0.326 




n (0a) 


0.575 




0.611 


0.611 






nm 


0.805 


0.920 


0.888 


0.888 


0.869 


0.858 


r>(ia) 


1.092 




1.166 


1.222 


1.141 






1.379 




I AAA 




1.467 














1.738 






1.667 


1.839 


1.721 


1.499 


2.010 


1.448 


Q(IIa) 


1.897 






1.777 






r.("b) 


2.184 












Q( m ) 


2.471 


2.759 


2.554 


2.110 


2.879 


1.985 


n (IIIa) 


2.759 






2.388 


3.205 




n (iiib) 


3.046 






2.721 






Q(IIIC) 








2.999 








3.276 


3.678 


3.443 


3.609 


3.748 


3.969 




4.138 


4.598 


4.331 


4.220 







